Targeted mixing in an array of alternating vortices 
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Transport and mixing properties of passive particles advected by an array of vortices are investi- 
gated. Starting from the integrable case, it is shown that a special class of perturbations allows one 
to preserve separatrices which act as effective transport barriers, while triggering chaotic advection. 
In this setting, mixing within the two dynamical barriers is enhanced while long range transport is 
prevented. A numerical analysis of mixing properties depending on parameter values is performed; 
regions for which optimal mixing is achieved are proposed. Robustness of the targeted mixing prop- 
erties regarding errors in the applied perturbation are considered, as well as slip/no-slip boundary 
conditions for the flow. 

PACS numbers: 05.45.Gg, 47.52. +j, 47.51. -|-a 



I. INTRODUCTION 

Since its uncovering, chaotic advection has drawn 
much attention as its consequences on both transport and 
mixing properties of a given flow are fundamental. This 
phenomenon is intimately related to Lagrangian chaos, 
and translates the fact that despite the laminar charac- 
ter of a flow Lagrangian trajectories of a fluid or pas- 
sive particles may end up being chaotic. In this setting, 
transport and mixing properties are drastically changed 
In chaotic regions of the flow, mixing induced 
by molecular diffusion becomes often negligible in regards 
to the mixing induced by the dynamics. Regarding trans- 
port properties, the triggering of chaotic advection also 
plays a key role. Indeed, in contrast to the fully predic- 
tive integrable situation, tackling transport properties of 
individual passive particles is subject to sensitivity to ini- 
tial conditions and implies a necessary probabilistic ap- 
proach. This leads in some cases to a diffusion equation 
with an enhanced diffusion coefficient when compared to 
the molecular diffusion one , but also to non-Gaussian 
properties of transport (see for instance [1, Q), as for in- 
stance super-diffusive transport. One then often resorts 
to modeling transport using a fractional diffusion equa- 
tion . All these properties have drawn much attention 

not only for its impact o n g eophysical flows and magnc- 
■ ■ ■ -r^r-^r-^^[f||g.r-.r-^r-^ - - 



tized plasmas @, [ll [Tl [Tl [Tl 111 [M [13, [111 , but also 
because of the potential applications of its enhanced mix- 
ing properties in chemical engineering and micro-fluidic 
devices 

Lagrangian chaos in two dimensional incompressible 
flows is triggered gcncrically when the flow becomes time- 
dependent [2, 121|,1^123|,12J|. The trajectories of fluid par- 
ticles or passive tracers are not confined on field lines and 
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chaos appears. For these type of flows the dynamics of 
fluid particles is Hamiltonian, with the stream function 
acting as the Hamiltonian. The canonically conjugate 
variables are the space variables, and the phase space 
is the two dimensional physical space. This particular- 
ity allows direct visualization of phase space in experi- 
ments, making it a test bed used to confront theoreti- 
cal results on Hamiltonian dynamics with experiments. 
More specifically, passive tracer dynamics belong to the 
class of Hamiltonian systems with one and a half degrees 
of freedom. General mixing and transport properties of 
these systems are now well understood, especially when 
the time dependence is periodic and Poincare maps are 
computed to analyze phase space topology. Typically 
the dynamics is not ergodic : A chaotic sea surrounds 
various islands of quasi-periodic dynamics. The anoma- 
lous transport properties and their multi-fractal nature 
are then linked to the existence of islands and the phe- 
nomenon of stickiness observed around them [25L |26| , 
while mixing is enhanced in the chaotic sea but has to 
rely on molecular diffusion in regular regions. 

In most studies regarding this type of phenomena, the 
time dependent perturbation is given a priori or self- 
generated. Transport and mixing properties are thor- 
oughly investigated and general laws are extrapolated or 
the origin of phenomena explained (see [26l. [271 [28l. l29j ) . 
The influence of phase space topology and its under- 
standing is clear, and can be used to explain synchro- 
nization phenomena [3(| . However one is still somewhat 
constrained by an a priori imposed time dependence. Re- 
cently, approaches of tailoring specific perturbations in 
order to modify phase space have been proposed [sil [3^ 
and a specific one has been applied to an array of al- 
ternating vortices [s^. Due to the strong influence of 
invariant tori forming regular islands on global transport 
and mixing properties, acting on phase space topology 
even locally (for instance, by building a transport barrier 
[31I [ssj or by destroying regular islands [33]) can have 
strong consequences. 

In this paper we address the problem of targeted mix- 
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ing, a work already started in Ref. (33|. We refer to 
targeted mixing as the process of leveraging only one of 
the consequences of chaotic advection, namely enhanc- 
ing mixing, while containing particle transport within a 
finite sub-domain of phase space. For this purpose we 
consider the dynamics of passive particles in an array 
of alternating vortices and tailor phase space in order 
to achieve the desired property. From the experimental 
point of view, acting on phase space is often not easy, be- 
cause one has to act also on particle momenta in general. 
However this may be less of a problem in the context 
of two-dimensional incompressible flows due to the du- 
ality between physical space and phase space. Moreover 
mixing within flows has a tremendous number of appli- 
cations. The primary interest in the flow of an array of 
alternating vortices resides in the fact that being gen- 
erated by quite a few hydrodynamic instabilities, it may 
be considered as one of the founding bricks of turbulence. 
This flow is easily accessible to experimentalists for in- 
stance using magnetohydrodynamics techniques similar 
to Rayleigh-Benard convection with a control over the 
flow [21I, [2^, [l^l ■ As such, understanding its influence on 
the advection of passive or active quantities is considered 
a necessary first step in order to uncover the different 
mechanisms governing transport in general or reaction- 
diffusion processes such as front propagation in turbu- 
lent flows JMj, j35„ ,36^ 38] . The stream function which 
models an experiment in a channel with slip boundary 
conditions writes: 



^'o (x, = sin X sin y, 



(1) 



where the x-direction is the horizontal one along the 
channel and the y-direction is the bounded vertical one. 
The dynamics given by the stream function ([T|) is inte- 
grable and passive particles follow the stream lines, no 
mixing occurs. In the experiments, a typical perturba- 
tion f{x,y,t) is introduced as a time dependent forcing 
in order to trigger chaotic advection and then to study 
the resulting transport and mixing properties. More pre- 
cisely we consider perturbations which modify the stream 
fimction as: 



*c(a;,y,t) = *o(a; + /,y). 



(2) 



We show how to identify the perturbations / which 
preserve transport barriers and at the same time, enhance 
mixing properties. 

The paper is organized as follows: In Sec. |TT] we re- 
call some basic notions of passive scalar dynamics in two 
dimensional incompressible flows. Then we give a short 
review of possible Hamiltonian control techniques we use 
in order to tailor a perturbation best suited for our needs. 
In Sec. mil we apply these techniques to flows modeled by 
the stream function ^ . First we derive the perturbation 
needed in order to build virtual barriers along the chan- 
nel and thus limit transport within only a small region 
of phase space. Then among all possible perturbations, 
we define criteria needed to achieve good mixing within 
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FIG. 1: Schematic representation of the channel with a trans- 
port barrier located at s = ipd/). 



the two barriers and analyze for which perturbations and 
parameters these criteria are satisfied. Then, we consider 
the robustness of the proposed perturbation with respect 
to a simpler time dependence, slip boundary conditions, 
three dimensional effects and molecular diffusivity. After 
showing the efficiency of such perturbations in enhanc- 
ing mixing, we analyze the set of parameters for which 
efficient mixing is expected using the residue method. 



II. A STRATEGY FOR MIXING INSIDE A 
CELL 

The term advection by definition, relates to the action 
of being moved by and with a flow. In mathematical 
terms, this translates in the general equation for a passive 
tracer: 



r = v(r, t) 



(3) 



where r locates the passive tracer in space, v corresponds 
to the velocity field, and r corresponds to the time deriva- 
tive of r. In the case of two dimensional incompressible 
fiows, the study of the field can be described by a scalar 
function, i.e the stream function ip{r,t). The velocity 
field is then obtained by v = curl(^ z), where z is the 
unit vector normal to the flow. Equation ([3]) is rewrit- 
ten using the stream function and exhibits a Hamiltonian 
structure for the flow : 



dy ' 



y 



9^- 

dx 



(4) 



where (x, y) corresponds to the coordinates of the tracer 
on the plane. The space variables {x,y) are canonically 
conjugate for the stream function ^ which acts as the 
Hamiltonian of the system. Hence the phase space is 
formally the two dimensional physical space (with the 
addition of time). 

In this section, our aim is to briefly recall some Hamil- 
tonian techniques allowing to some extent, to tailor phase 
space by adjusting appropriately the perturbation / as 
in Eq. ^ and its parameters. These techniques will be 
subsequently used to achieve targeted mixing for the con- 
sidered flows, namely construct a perturbation in order 
to obtain optimal chaotic mixing in a localized region of 
phase space. First we construct a family of perturbations 
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which create transport barriers, then we use the residue 
method [s^ in order to find optimal mixing regimes in 
parameter space after the barriers have been created. 

A. Constructing a perturbation with a barrier 



In order to have an invariant curve acting as a barrier in 
the channel, the perturbation / is such that the stream 
function evaluated at a; = ipijj) is only a function of 
time, since the barrier is time- independent (in the moving 
frame). For simplicity, we consider solutions for which 
the stream function vanishes : 



We consider a generic time-independent stream func- 
tion (Hamiltonian) 5*0(2;, y) describing a fluid in a two- 
dimensional channel of height tt, i.e. {x, y) € Mx [0, it] (see 
Fig. [T]) . We assume that there exists an invariant curve 
which prevents the advection of tracers from the left to 
the right of that curve. We denote the equation of this 
curve X = (p{y). The invariance condition [x = (p'{y)y if 
X = (p{y)] of tracers on this curve translates into : 

^0 {f{v),y) = constant, 

for all y £ [0, tt] by using Eqs. dH). In the case of Eq. ([1]), 
this transport barrier is a heteroclinic connection be- 
tween two hyperbolic fixed points located at y = and 
y ^TT. 

In order to enhance mixing, the system is perturbed by 
a time-dependent forcing (depending also on x and y in 
general). As a consequence, its dynamics is generically no 
longer integrable, and chaotic trajectories fills portions 
of the channel (these are parts where mixing occurs). 
However, as a side effect, the transport barrier is broken 
and trajectories start to diffuse along the channel (in the 
x-direction) . 

In the following, we propose to design a time- 
dependent forcing of the fiow described by which pre- 
serves the transport barrier as well as the chaotic mix- 
ing. The aim is to find an appropriate forcing which pre- 
serves a bounded domain of the channel with an enhanced 
chaotic mixing inside. This domain can be bounded by 
two dynamical barriers as mentioned above. 

A main practical requirement is that the forcing should 
be as simple as possible to be implemented. Conse- 
quently, we start by investigating perturbations which 
only depends on y and t. Furthermore, for practical rea- 
sons, we restrict the perturbations to time-periodic ones, 
and the period is chosen as 2tt without loss of generality. 
More precisely, we search for a perturbation f{y, t) which 
modifies the stream function into : 



■^c{x,y,t) = *o [x + J{y,t),y]. 



(5) 



In order to simplify the computations, we perform a 
generic translation in x (which corresponds to a canonical 
transformation in the Hamiltonian setting) : 

X = X + dyf3{y,t), 

y = y, 

where [3 is & function to be specified later. In the new 
variables x and the dynamics is described by the 
stream function : 

4'e(i, y, t) = *o [i + /(y, t) - dyp{y, t),y]- dtp{y, t). 



*o My) + /(y, t) - dyPiy, t),y]- dtf3{y, t) = 0, (6) 

which is a single equation with two unknown functions, 
f3 and /. In principle, there are an infinite set of solu- 
tions. Depending on other requirements, some solutions 
are more appropriate than others. In the following, we 
choose 

fiy,t) + ^iy)~dyPiy,t)^^t), 
where $(t) is any function of t. Equation ([6]) implies 

dtp{y,t)^^o [<&(*), y]. 

This equation has solutions provided that the mean value 
of "ifo [^{t), y] with respect to time vanishes for all y. This 
guides the choice for the function $. A possible solution 
for f3 is 



/3(y,t) = r*o [$(i),y] 



where 



fe#0 

for V = J2k '^k{y)&^'^* ■ The perturbation / is given by 

/(y, t) = m - ip{y) + Tdy^^ m), y] ■ (7) 

The stream function ^'c given by Eq. ^ has an invariant 
curve whose equation is 



x = ^{y)~Tdy^^mt),y). 



(8) 



There exist many more solutions / than the one indi- 
cated here. These solutions are obtained by using more 
complex translation functions (3. However, they contain 
a lot more Fourier modes in y, which make them more 
difficult to implement for the cases we consider in Sec. lIIII 
We notice that the perturbation / given by Eq. d?]) as well 
as the equation of the transport barrier ([8]) are parame- 
terized by a time-periodic function <I>(i). This function $ 
is in general parameterized by essentially two parameters, 
its frequency and its amplitude. In parameter space, the 
dynamics characterized by the stream function ^'c exhibit 
drastically different behaviors depending on the values of 
parameters. In what follows, we use periodic orbits to de- 
termine the regions in parameter space where complete 
mixing inside the transport barriers occurs. This method 
is briefly explained in the next section. 
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B. Periodic orbit analysis 

The dynamics of the perturbed flow as the one given 
by the stream function ^'^ can be investigated by looking 
at the (linear) stability of specific periodic orbits, using 
indicators such as the residue. In order to "control" the 
flow, one can monitor the residues by varying the pa- 
rameters of the system (like the frequency and the am- 
plitude of the forcing), until specific bifurcations occur. 
In particular, the residue method allows one to predict 
the break-up (or creation) of invariant tori, which entails 
an enhancement (or reduction) of chaotic mixing. 

We consider a Hamiltonian flow with 1.5 degrees of 
freedom which depends on a set of parameters A G M'" : 

r = JVr*(r,i;A), 

where r = (x, y) and J = ^ ) ' "'"^ order to analyze 

the linear stability properties of the associated periodic 
orbits, we also consider the Jacobian J*(r) which evolves 
according to the tangent flow written as ^] : 

lj*(r)=Jv2vI/(r,i;A)J*, (9) 

where is the two-dimensional identity matrix and V^^ 
is the Hessian matrix (composed of second derivatives of 
\1/ with respect to its canonical variables). For a given 
periodic orbit with period T, the linear stability prop- 
erties are given by the spectrum of the two-dimensional 
monodromy matrix j'^ . These properties can be syn- 
thetically captured in Greene's definition [13, [Hj of the 
residue : 



since det j'^ ~ 1. In particular, if R e]0, 1[, the periodic 
orbit is elliptic; if i? < or i? > 1 it is hyperbolic; and if 
R — and i? = 1, it is parabolic. 

The residues of a set of well selected periodic orbits 
provide - through linear stability analysis - information 
to detect the enhancement as well as the reduction of 
chaotic mixing [H, \^ . The actual change of dynamics 
has to be checked a posteriori by a nonlinear stability 
analysis (a Poincare section for instance). The residues 
are used to discard regions in parameter space where 
large elliptic islands are present. 

For this purpose, an alternative strategy is to use a 
"brute force" method and scan a whole range of physi- 
cally relevant parameters, analyze transport and mixing 
properties (by a Poincare map inspection for instance) 
and conclude on a domain of parameters where optimal 
mixing is achieved. Though it would be a complete anal- 
ysis, this strategy is not reasonable to adopt because of 
the high number of cases to consider and also the com- 
puter time it takes to analyze a single case. 

In order to circumvent these difficulties we choose to 
consider the residue method described in this section and 



follow in parameter space the stability of a well selected 
set of periodic orbits. Indeed due to the Hamiltonian 
nature of passive particles, one can expect a direct corre- 
spondence between non-mixing regions in physical space 
and islands of stability in phase space. The higher the pe- 
riod of the island the smaller is its size, hence by following 
the linear stability of elliptic periodic orbits with small 
period in parameter space, one should be able to define 
a potential optimal mixing region for which these orbits 
are unstable and associated with a mixing enhancement. 
Once this set of main periodic orbits has been defined by 
close inspection of several situations, the mixing region 
is obtained by looking at the bifurcation curves in pa- 
rameter space, e.g., the set of parameters such that the 
residue is equal to one [3^. We will show in the next 
section how to combine the creation of transport barriers 
and the residue method in a particular example of stream 
function given by Eq. ([T]). 

III. ACHIEVING TARGETED MIXING IN AN 
ARRAY OF VORTICES 

The stream function given by Eq. ^ models a cellular 
flow consisting of alternating vortices. If we restrict our- 
selves to y e [0, tt], we have a channel of alternating vor- 
tices with slip boundary conditions. From the Hamilto- 
nian perspective, the advection of passive tracers is given 
by Eq. ([4]). Since the flow is steady, trajectories coincide 
with the fluid streamlines depicted in Fig.[2Ja). 

Few basic facts explain the structures of the dynam- 
ics : Boundary conditions given in y = and y = tt 
constitute invariant curves. The system is 27r-periodic in 
the x-direction, i.e. along the channel. The system has 
hyperbolic flxed points located at x = rmr for m G Z and 
y ~ or y ~ n. These points are joined by vertical het- 
eroclinic connections for which the stable and unstable 
manifolds coincide corresponding to roll boundaries, at 
the origin of the cellular structure of the flow. 

In order to obtain chaotic advection for two dimen- 
sional flows, one has to perturb the flow by a time- 
dependent forcing. For example, one can periodi cally 
force the roll patterns to oscillate in the cc-direction [2l| , 
in which case the stream function reads : 

\I'i(x, y, i) — sin(2; + e sinoji) sin y, (10) 

where f{x, y,t) = e ainujt acts as a perturbation. The pa- 
rameter e and uj are respectively the amplitude and the 
angular frequency of these lateral oscillations. The barri- 
ers (heteroclinic connections between the hyperbolic pe- 
riodic orbits) are broken under the perturbation /, which 
cause the passive particles to undergo chaotic advection 
along the channel [i^ . 

The streamlines corresponding to the stream function 
^I'l given by Eq. PU)) are depicted in Fig. ^a). We ob- 
serve that the structures are essentially the same as those 
of the stream function ^'o given by Eq. ([1]). However, the 
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periodic forcing now drives back and forth the roll pat- 
terns in the x direction with a period ^-n ju. The surfaces 
y — Q and y — are left invariant by the perturba- 
tion and the hyperbolic orbits persist on these surfaces. 
The Poincare section for w = 1 and e = 0.8 (see Fig. [2] 




X 



FIG. 2: (a) Streamlines at t = and and (6) Poincare section 
of the stream function (|10p . The parameters are lo = \ and 
e = 0.8. 

(6)) shows how passive particles are spreading along the 
channel. Particle transport from roll to roll is greatly 
enhanced. However, an unmixed area characterized by 
regular trajectories is still present at the center of each 
vortex : It is composed of invariant tori of the dynam- 
ics. Consequently, though diffusion has appeared in the 
system due to the time-dependent perturbation, some 
regular patterns are persistent. Note that increasing the 
amplitude e does not make them disappear; in particular, 
in the limit of large e (with lo fixed) , the system will be 
integrable as well. 



A. Building transport barriers 

In what follows, we propose to adjust the time peri- 
odic forcing as explained in Sec. Ill Al First we notice that 
^o(™7i'7 y) — ^ for any y £ [0, tt] and all m e Z. We ap- 
ply the construction of the perturbation with ^(y) — 
(or equivalently (/5(j/) = 2m7r). In Eq. ([7]), we choose 
<I>(t) = esinwi which gives the perturbation : 

/(y, i) — esinwi -I- lo~^ cosyCe{ujt) , (11) 

where 

C^{t) = rsin(esini) 

= -2 V:^^J72„+i(e)cos(2n-Hl)i, (12) 

n>0 

and J'n (for n G N) are Bessel functions of the first kind. 
In the numerics, we truncate the sum p2p to five modes. 
Very similar results are obtained with a higher number 
of modes. 

We notice that since the stream function is still 2tt 
periodic in the x-direction by construction, the invari- 
ant surface which has been created around a; = is also 



translated around 2rmr, m G Z. The equations of these 
transport barriers along the x-direction are : 

X ~ 21717: — uj~^ cosy Ce{ujt), (13) 

for all TO G Z. Notice that these barriers exist for ar- 
bitrary values of lo and e. Each of these barriers are 
heteroclinic connections between two hyperbolic periodic 
orbits: 

x{t) — 2rmT — uj^^Ceiujt) at y = 0, 
x{t) — 2rmT + uj^^Ce{ujt) at y = tt, 

which move in opposite directions. Furthermore, the top 
and bottom boundaries of the channel remain invariant. 
This comes from the fact that the perturbation given by 
Eq. (fTT|) is only applied in the x term of ^'q. 

Figure[3Ja) depicts the streamlines of the stream func- 
tion 

\['c(a;, y, t) = sm[x + es'mojt + uj~^ cosyCf{ujt)] s'my, 

(14) 

a,t t — for u — 1 and e ~ 0.8. We notice that the 
displacement of the rolls remains periodic and parallel to 
the x-direction with an additional oscillating shear. The 
Poincare section depicted on Fig.[3{6) reveals the dynam- 
ics of tracers which is very different from the one given by 
the stream function given by Eq. (fTO)) (see Fig. [2^6)). 
The main differences are that barriers suppressing long 
range chaotic transport along the channel are restored 
around a; = (mod 27r) (bold curves), and that efficient 
mixing is achieved within the cell confined by two barri- 
ers. We observe that passive particles appear to invade 
the whole confined cell until the fluid is apparently fully 
mixed; most of the regular trajectories observed for the 
stream function ^'i are broken by the perturbation. 



(a) ^ , (b) 




FIG. 3: (a) Streamlines at t = and (b) Poincare section 
of the stream function (|14p . The parameters are ui — 1 and 
e = 0.8. 

In Fig. |4l a numerical simulation of the dynamics of a 
dye of tracers in the fluid is shown. The left column de- 
picts the dynamics of tracers for the stream function (jlOp . 
The right column shows the mixing of a dye within a cell 
delimited by two barriers created by the stream func- 
tion (fH)) . We see that the scattering of the dye, which 
leads to mixing, occurs through a combination of stretch- 
ing and folding of the dye in both cases. 



FIG. 4: Numerical simulation of the dynamics of a dye at 
t = 4T, t = 6T, t = 18T and T = 2-k (from top to bottom) : 
left column for the stream function pup and right column for 
the stream function (|14p . The parameters are uj = 1 and 
e = 0.8. 



B. Mixing analysis: Local Lyapunov exponent and 
mean recurrence time analysis 

As the absence of the regular islands from the Poincare 
sections does not guarantee the homogeneity of mixing, 
the study of local properties of mixing in phase space 
may provide extensive information. For this purpose we 
consider two different types of analysis, namely the Lya- 
punov map and the mean recurrence time. Both analysis 
are performed within the space of initial conditions. 

First, in order to get insight into the action of the per- 
turbation on the local stability properties of the system, 
we compute the Lyapunov map. This method provides 
local information in phase space. It has been introduced 
to detect ordered and chaotic trajectories in the set of 
initial conditions. It associates a finite-time Lyapunov 
exponent u with an initial condition (xojj/o) i — 0. 
Let us consider the tangent flow ^ , and define the max- 
imum finite-time Lyapunov exponent by integrating the 
flow and the tangent flow over some time r starting with 
some initial condition (xq , yo ) : 

i/(xo,?;o,t) = -log|Amax(a;o,2;o,T)|, 

T 

where Amax(2;o, J/O: ''') is the largest (in norm) eigenvalue 
of J"^ (one can also use the eigenvalue of J'^*J'^ where 
J'^* is the transposed matrix of J"^). 

From the inspection of the map {xo,yo) ^ v(xo,ya,T) 
for some given time r, one distinguishes the set of ini- 
tial conditions leading to regular motion associated with 
a small finite-time Lyapunov exponent, from the chaotic 
ones with larger finite-time Lyapunov exponents. Hence, 
this map reveals the phase space structures where the 



FIG. 5: Lyapunov maps for the dynamics of tracers given by 
the stream function given by Eq. (I14p for cu — 1.67 and 
e = 0.63 (upper panel) and for ui = 1 and e = 0.8 (lower 
panel) at a time T — 2007r. 

motion of tracers is trapped on invariant tori, i.e. they 
highlight islands of stability located around elliptic pe- 
riodic orbits. Mixing regions are characterized by high 
values of their finite-time Lyapunov exponents. 

Figure [5] represents the Lyapunov maps for the dy- 
namics of tracers given by the stream function (fT4|) for 
UJ = 1.67 and e = 0.63 (upper panel) and for lu — 1 
and e — 0.8 at a time T = 2007r. The dark regions are 
characteristic of small values of the Lyapunov exponent. 
We notice that Fig. [5] (upper panel) shows small remain- 
ing islands which are barely noticeable in the Poincare 
section (see Fig. 3 of Ref. [33[). For mixing studies, the 
Lyapunov diagnostic seems to be an appropriate tool to 
reveal small non-mixing regions. These regular regions 
have disappeared for w = 1 and e = 0.8 (lower panel). 
However, the Lyapunov map is not able to identify the 
transport barriers which means that locally near the bar- 
riers, the motion is as chaotic as inside the cell. 

As seen above, by looking at the Lyapunov map, one 
can infer local mixing properties of the fiow. However 
one can notice that since the created barrier is a separa- 
trix and not a KAM torus as for instance in Ref. [31| , the 
existence of the barrier cannot be detected by the Lya- 
punov map. To complement this analysis, we consider a 
second diagnostic namely a recurrence time analysis. An 
interesting property of return time distributions stems 
from the fact that they are known to be sensitive both 
to local and global dynamical properties of phase space. 
For instance, being in the neighborhood of a hyperbolic 
periodic orbit versus an elliptic one should affect the dis- 
tribution Therefore the distribution should be af- 
fected if computed in the neighborhood of a separatrix, 
or if trapped within a regular region. Since we have al- 
ready analyzed the local stability properties of the flow by 
computing the Lyapunov map, we will only consider the 
average first return time and define a scalar field r(x,?/) 
in the space of initial conditions. One can indeed assume 
that if near a given point (xo^yo) with positive i^(a;o,?/o) 
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FIG. 6: Mean recurrence time in initial condition space of 
tracers given by the stream function given by Eq. (|14|) 
for uj — 1.67 and e = 0.63 (upper panel) and for ui = 1 and 
e = 0.8 (lower panel). Initial conditions space is divided into 
100 X 50 cells. Trajectories are computed for ~ 10® periods, 
time steps is dt = T/200. . 

the average return time T(xo,yo) is large, then a dye of 
fluid would explore a large part of phase space and so it 
would be best to drop the dye in its neighborhood, than 
in an other point with similar v but a smaller t. 

Typically when analyzing return time statistics from a 
numerical perspective, one defines a small area T in phase 
space and compute a distribution of return times of tra- 
jectories leaving the area. On one side due to the fact 
that, for Hamiltonian systems with finite phase space, 
the average return time is finite and scales as l/T (Kac's 
lemma), one has to be cautious not to take F too small 
in order to carry long enough simulations and capture 
enough events to build a characteristic distribution. On 
the other hand, return time distributions are supposed to 
be computed for T ^ 0. For the considered flow, due to 
symmetries we consider [0, tt] x [0, tt] as the space of ini- 
tial conditions which has been divided regularly in 2500 
small squares. For each square, the mean return time has 
been computed using two trajectories computed for 10^ 
periods. Given the size of F, we collect for each cell about 
300 events. The results are presented in logarithmic scale 
in Fig. [SI The parameters have been chosen identical to 
the ones used in Fig. [SI In contrast to the Lyapunov map, 
one sees that this diagnostics finds the barriers (region 
with long return times). For lo ~ 1.67 and e — 0.63 (up- 
per panel), one can also see small regions with quite low 
return times : They correspond to small regular islands 
as mentioned previously. The small return time regions 
have disappeared for to — 1 and e — 0.8 (lower panel). 

Hence, in order to get an accurate picture of the mixing 
properties of the cell, one has to combine the information 
of both the local Lyapunov exponent and the local av- 
erage return time, for example by computing the scalar 
field vItojI/o) X T(xo,yo)- Indeed this map shall give us 
information on good mixing regions, and provide as well 



the location of where to drop initially the dye to achieve 
a faster homogenization and mixing. 

C. Robustness 

We have identified a family of perturbations which, 
while keeping the cellular structure of the flow, enhance 
considerably mixing properties within the cells. The per- 
turbations have been constructed for a very specific flow; 
hence a natural question arises in practical situations: 
whether or not the considered perturbed system is robust 
with respect to small changes or errors in the applied per- 
turbation. Indeed robustness is a key concern in order for 
an experimental set up using this type of perturbations. 
Below we analyze robustness with respect to four factors: 
the truncation of the time series giving the time depen- 
dence of the perturbation, the slip boundary conditions, 
three dimensional effects, and molecular diffusivity. 

1. Truncation of the series 

From the experimental perspective we may expect 
some difficulties in implementing the whole series C^(t) 
given by Eq. (jl2p . One may wonder how the barrier 
and mixing properties are affected when one truncates 
the series and retain for instance only the first term of 
the perturbation term (|12p . meaning only one mode is 
used for the perturbation, Ce is replaced by —2Ji (e) cos t. 
Figure [7] represents the plot of these two functions for 
uj = 1 and e = 0.8. The small discrepancy between both 
functions might affect significantly the transport barriers 
since it is well known that heteroclinic orbits are very 
sensitive to perturbations and are generically broken by 
an arbitrarily small perturbation. 

Trajectories of passive particles with a dynamics given 
by the stream function 

^l°'\x,y) ~ sin[a;+esincjt — 2cj^"^Ji(e) cos cosy] siny, 

(15) 

are displayed in Fig. [5] One can see that the barrier 
is leaking while mixing properties do not seem to be af- 
fected significantly. One has to notice that given the time 
length and the amount of passive particles considered, 
the leak is small. The truncated series is still achieving 
a good "targeted mixing" . 

2. No slip boundary conditions 

If one wants to take into account the thin bound- 
ary layers present on the boundaries of the channel, it 
amounts to consider that the fluid is confined between 
the two surfaces y = and y = it with no-slip boundary 
conditions. In this case, the stream function is modified 
into : 

■^o{x,y) = sinx W{y) , (16) 
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FIG. 7: Functions C^{t) (solid line) and -2Ji(e) cost (dotted 
line) for t e [0, 27r] and e = 0.8. 




FIG. 9: Poincare section of the perturbed stream func- 
tion (|18p . The parameters are = 1 and e = 0.8. The 
integration time is 1000 periods. 
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FIG. 8: Poincare section computed with the trajectories of 
1000 passive tracers and dynamics given by the stream func- 
tion (|15p . The parameters are a; = 1, e = 0.8. The integration 
time is 1000 periods. 

where no-slip boundary conditions x — y — Oa,ty = 
and y = TT are obtained with 

W{y) = cos{qoy) - Ai cosh{qiy) cos{q2y) 

+ A2sinh{qiy)sin{q2y), (17) 

with y ^ y/n ~ 1/2, qo = 3.973639, qi = 5.195214, 
q2 = 2.126096, Ai = 0.06151664 and A2 = 0.103887 
(see Ref. [45|). In this setting the resulting streamlines 
of the unperturbed stream function (jl6p are very simi- 
lar to the ones in Fig. [D^a) with the difference that the 
velocity vanishes at the top and bottom of the channel. 

In order to test the effect of the proposed perturba- 
tion (|lip on the stream function (fT6|) . we consider the 
following perturbed stream function : 

^l{x, y, t) — sin[x -|- esmut + uj'^^ cos j/Ce(u;t)] W{y) . 

(18) 

Note that the no-slip boundary conditions are maintained 
at 2/ = and y = tt. A poincare section for the stream 
function is displayed in Fig. [H] for cj = 1 and e = 
0.8. The barriers are broken and we observe long range 
chaotic transport of passive particles along the channel. 
However, the mixing properties are maintained in the 
major part of the channel, except for some areas where 
some very small regular regions remain. 

Nevertheless, the exact perturbation of the stream 
function ^TE\i can be derived following the method devel- 



FIG. 10: Poincare section of the stream function (|19|) . The 
parameters are u = 1 and e = 0.8. 

oped in Sect. Ill Al It gives the following stream function : 

*c(a;, y, t) = sin[a; + esmujt + uj^^W' {y)Ce{ujt)]W{y), 

(19) 

which preserves the no-slip boundary conditions at y = 
and y = TT. The resulting Poincare section depicted 
on Fig. [ini for a; = 1 and e — 0.8, shows, as in the 
case of Fig. EJb) , that the stream function keeps the 
transport barriers and the mixing properties. 

3. Three dimensional effects 

If the flow is bounded, depending on the size of the 
boundary layers and the fluid extension in the third di- 
rection, non uniform vorticity can give rise to a secondary 
instability leading to a three-dimensional flow (Ekman 
pumping). Hence, for the considered perturbed flow, we 
may have to take into account the weakly three dimen- 
sional case, which may be given by the empirical flow 
(see for instance [46j): 

X = — sin a; cosy + £1 sin 2a; sin z, 

y — cosxsiny -|- ei sin2?/sinz, 

z ~ 2ei cos z[cos 2.x + cos 2?;]. (20) 

Note that the strength of the third component of the flow 
is characterized by ei. 
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In order to study how well the two-dimensional bar- 
rier fares in this three-dimensional flow, we apply the 
perturbation given by Eq. ([TT|) to the right hand side of 
Eq. (pni) . The perturbed system is then given by: 

i = — sin cosy -|- ei sin sin z 

+Lj~^Ce{ujt) COS Xs sin^ y, 
y = cosa;s siny -|- ei sin2ysin2;, 
z ~ 2ei cos z[cos 2a;s -)- cos 2y], (21) 

whereas = x + esmut + uj^^ cosyCe{ujt). We notice that 
an additional term has been added to x in order to ensure 
a divergence free field. 

Keeping a two-dimensional point of view of the sys- 
tem (HU, we visualize the projections in the {x,y) plane 
of the position of passive tracers. When considering 
ei — 0.005 and 4.10'' passive particles at time t = lOT, 
where T = 2n, an effective barrier remains as it is shown 



in Fig. 11(a) and mixing properties are not affected, but 



we observe some advected particles which escape from 



the cell. However as depicted in Fig. 11(b) when the in- 
tegration time is t = lOOT (for the same value of ei and 
number of particles), the barrier still influences the mo- 
tion but leaks since a more significant number of particles 
get through these barriers. 




(a)t = lOT 



(h)t = WOT 



FIG. 11: Projection of trajectories of the system (|2U|I for lu = 
1, e = 0.8 and ei = 0.005. 
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FIG. 12: Trajectories of 1000 passive tracers given by 
Eq. (|22|l . (top). Local density of tracers (bottom). The 
parameters are fj, — 4.10"^ (Peclet number Pg ~ 600), lo — 1 
and e = 0.8. The integration time is 1000 periods. 



where /i = 4.10"^ (which corresponds to Peclet number 
of Pe « 600), oj = 1 and e ~ 0.8. In order to avoid cross- 
ing across the "walls" y = tt and y = 0, we took by = 0. 

One can see that for this type of Peclet values, the 
transport barrier is indeed leaking, and in fact the trun- 
cation of Ce(t) does not affect how much the barrier leaks. 

In summary the robustness of the proposed pertur- 
bation in different settings has been investigated. One 
may infer that the most drastic effects are induced by 
boundary conditions. On the other hand provided that 
the right perturbation is computed, one concludes that if 
three-dimensional effects are weak and Peclet number is 
high enough, a truncation to the first term of the series 
Ce(t) is sufficient, and actually implementing more terms 
seems useless. 



D. Targeted mixing regimes 



4- Molecular diffusivity 

Finally, we have up to now considered ideal passive 
tracers, which are not subject to any molecular diffusiv- 
ity, which may be a good approximation for high Peclet 
numbers. However it is likely that for any finite molecular 
diffusivity the barrier will leak. In order to illustrate this 
phenomenon, we consider that tracers are actually sub- 
ject to a Langevin equation associated with the stream 
function c given by Eq. ^4]) : 

X = -^+b,{t), y^—l+by{t), (22) 

where b,j.{t) and by{t) are two independent delta- 
correlated white noises, with zero mean and a given am- 
plitude fi. Numerical results are displayed in Fig. 1121 



It was shown in Sec. Ill Al that a perturbation allows 
one to localize tracers into a finite volume of phase space. 
However, depending on the values of the parameters of 
the perturbation (w, e) in Eq. pi)) , some regular islands 
may exist and prevent complete mixing inside the cell. 
In this section we propose to identify the domain of pa- 
rameters (if any) such that these islands do not exist and 
the cell become fully mixing. 

This domain in parameter space can be determined by 
analyzing the linear stability of a few periodic orbits of 
the system as discussed in Sec. Ill Bl These orbits are 
those with low rotation numbers, around which the res- 
onant islands organize. Indeed, an island is organized 
around a central elliptic orbit : If the latter were to turn 
hyperbolic, the dynamics might become (locally) chaotic, 
since generically chaos is expected in the neighborhood 
of hyperbolic orbits by an infinite number of intersections 
between its stable and unstable manifolds. In order to 



10 



monitor these orbits, we use a scalar indicator of linear 
stability, such as the residue (see Sec. IIIBI) . The fully 
mixing regime will then correspond to the residues of 
selected periodic orbits being below or above 1. How- 
ever, the change of linear stability of a given periodic 
orbit might not change its nonlinear stability and still 
preserve invariant tori in its neighborhood. This linear 
stability analysis of periodic orbits has to be completed 
by an a posteriori check to determine whether or not the 
regular island has disappeared with the elliptic periodic 
orbit turning hyperbolic. 

Let us consider the following range of parameters 
{uj, e) G [0.6, 2.2] X [0, 2.2]. Inspection of several Poincare 
sections reveals that in this range, it is essentially the 
nature of eight periodic orbits which drives the mixing 
properties inside the cell. However, thanks to the sym- 
metry with respect to the point {x,y) — (7r,7r/2), this 
set reduces to only four orbits : Three of them have ro- 
tation number one, namely Oi, Of and Of, while the 
fourth one, called has a rotation number Q — 2. De- 
pending on the value of the parameters (w,e), these or- 
bits can be elliptic - and create resonant islands around 
them -, or hyperbohc. In order to describe the nature of 
each of these four orbits for a given value of parameters, 
we use the nomenclature [N{0i)N{0t)N{01)N{0^)], 
with N{0) = /i if O is hyperbolic, e if elliptic, and if 
it does not exist. Let us precise that there also exists 
a hyperbolic orbit with Q = 2, called O2 , which forms 
a Birkhoff pair with O^. As it remains hyperbolic in 
the range of parameters under consideration, it only pro- 
vides a better understanding of the bifurcation process, 
but does not influence the mixing properties inside the 
ceh. 

Figure fT3l depicts Poincare sections for eight different 
values of parameters. These cases illustrate some pos- 
sible mixing regimes in the considered range of parame- 
ters. Figures 13(a)p"3(d) (corresponding respectively to 
[hhOO], [heOO], [e/iOO] and [eeOO]) show how d and 
can coexist in their two forms (elliptic and hyperbolic), 
the full mixing regime being reached when both are hy- 
perbolic as in Fig. 13(a) where it is [hhOO]. We notice that 
when Oi or is elliptic, the mixing is only prevented 
by the elliptic island. Only small secondary islands are 
observed (for instance, in Fig. 13(c) ). This reinforces the 
importance of the considered set of periodic orbits for 
mixing properties. 

For smaller values of w, the orbits O" and Of may 
also appear (if w < 1.2 and large e for O", or w < 1 and 
small e for Of), and their nature have to be taken into 
account, as one can see on Fig. 13(e)| \hheO], where the 



hyperbolicity of Oi and is not sufficient to ensure the 
full mixing inside the cell, as Of is present in its elliptic 
form. Mixing can still be obtained in presence of Of , as 
can be seen on Fig. |13(f)| which corresponds to th e case 
[hhOh]. The opposite case [eeOe] is depicted on Fig. 13(h) 



where all three orbits are elliptic and almost no mixing 
occurs. 



Finally, one can see on Fig. 13(g) which is [eOOO] how 
for large values of w (typically beyond 2), no more 
exists, while Oi stays elliptic : Furthermore, new reso- 
nant islands have appeared, associated with new periodic 
orbits (with Q ~ 3 and Q = 5 as shown). This illustrates 
the fact that in another range of parameters, the dynam- 
ics may be guided by higher order periodic orbits which 
are associated with smaller islands. 



(a)aj = l,e = O.i 



{h)ui = 1.8, e = 0.5 



(c)u) = 1.2, e = 0.2 



{d)Lj = 1.6, e = 0.2 





• 

1 


I 
{c)uj -- 


2 4 6 8 
X 

= 0.58, e = 1.25 



1 



2 4 6 
X 

{f)uj = 0.58, e = 0.6 



ig)uj = 2.1, £ = 0.25 




(h)aj = 0.63, e = 0.01 



FIG. 13: Poincare sections in the different domains of param- 
eters {uj, e) for the stream function (|14l) . The stars, full circles, 
squares and triangles indicate respectively the locations of Oi , 
Of, of and Ot- 



A better insight in the bifurcation scheme can be 
gained by varying only one parameter at a time. Fig- 
ure [14] depicts the residue curves associated with the or- 
bits 01,0^,0^, Of and Of when varying e and keeping 
UJ constant. Two curves have been computed in the in- 
termediate frequency range co — 1.42 and w = 1.67. In 
the intermediate regime, only three of these five orbits 
exist, namely Oi (plain line), and (resp. upper 



and lower dash-dotted lines). 
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For w = 1.42 (see Fig. |14(a)| and for low e, both Oi 
and 0+ are elliptic : only partial mixing occurs. While 
O2 turns hyperbolic as e is increased (around e ~ 0.2), 
Oi remains elliptic until e « 0.4. Oi turns from ellip- 
tic to hyperbolic by merging with O^, according to the 
process : 



0^[h]+Oi[e] Oi[h]. 



(23) 



This bifurcation does not happen just for this par- 
ticular choice of parameters, but appears to be fairly 
generic in this range of parameters. For a better illus- 
tration of this process, the behavior of the eigenvalues 
of the monodromy matrix associated with Oi is shown 
on Fig. [15] (for cu — 1.67, similar to the case oj = 1.42). 
The two eigenvalues, initially conjugated on the unity 
circle (i.e. Oi is elliptic), encounter a period doubling 
bifurcation when they reach 1, and leave the unity cir- 
cle : Oi has turned hyperbolic. Then they are of the 
form (A, 1/A), A G M*. Eventually, the phenomenon will 
revert, the eigenvalues going back on the circle and Oi 
to ellipticity. For e between 0.4 and 0.9, the orbits Oi, 
O2 and O2 are all hyperbolic, acknowledging a com- 
plete mixing. Then, increasing further e, the orbit 
turns elliptic after e w 0.9, and hence mixing decreases. 
It remains so until it merges with the hyperbolic Oi, at 
e « 1.1, to give an elliptic Oi (parabolic at the transi- 
tion), according to the scheme : 



Ot[e] + Oi[h]^Oi[e]. 



(24) 



Beyond e w 1.4, the only remaining orbit, Oi, stays 
elliptic. 

For u = 1.67 (see Fig. 14(b) I, the dynamics is more 
regular : While bifurcation (|23p still occurs, resulting 
in an hyperbolic Oi, the residue of O2 never crosses 1, 
which means that the orbit stays elliptic. Around e ~ 0.9, 
it merges with the hyperbolic Oi to give an elliptic Oi, 
according to the bifurcation (|24p . Beyond this point, Oi 
will stay elliptic : Full mixing cannot be achieved for such 
values of uj. 

Three other curves illustrate the behaviors at low fre- 
quency, i.e. UJ = 0.8 and uj = 0.58, and at a high fre- 
quency UJ = 2.08. 

For UJ = 0.8 (see Fig. 14(c) I, the bifurcation scheme 



is more complicated due to the presence of the orbit Of 
and of (resp. right and left dashed line). The orbits can 
either exist in an elliptic or hyperbolic way, or not exist 
at all. The two latter cases, combined with Oi and O^ 
hyperbolicity, are suitable for a fully mixing regime. 

For low e, though O^ turns hyperbolic very soon (at 
e « 0.05), of is elliptic until e « 0.2, when it disappears. 
However, because of the ellipticity of Oi, full mixing can- 
not occur until it also turns hyperbolic, at e ^ 0.3. Then 
complete mixing is achieved until e « 1.5, when O^ turns 
back elliptic, soon followed by Oi when they merge ac- 
cording to the bifurcation scheme (IM)) . Note that around 



e « 1.6, orbit O" appears, but full mixing is no more pos- 
sible because of Oi ellipticity. 

Then, for uj = 0.58, the bifurcation scheme is more 
intricate due to the fact that Of turns hyperbolic. For 
low e, despite O^ and Oi soon turn hyperbolic (at e « 0.2 
for the latter), Of remains elliptic until e « 0.5, when 
its residues goes above 1 : full mixing is then achieved 
(see Fig. |13(f)[ ), until the orbit returns to ellipticity at 
e w 0.65. However, it disappears at e ~ 0.7, and then 
only Oi and O^ are present, in their hyperbolic form; 
thus mixing is achieved anew. It remains so until e w 1.1, 
when Of comes into play, being elliptic. Furthermore, it 
will also be hyperbolic (around e « 2), but only when Oi 
has turned back to ellipticity : mixing will not occur any 
longer. 



For large uj (uj = 2.08, see Fig. 14(e)), the dynamics 
is more regular. The orbits O^ and O2 do not exist, 
and the only remaining orbit, Oi, never encounters any 
bifurcation, and stays elliptic. Complete mixing cannot 
be achieved in this case either. 

Now, instead of varying the amplitude e of the forcing, 
we vary its frequency uj and keep e constant. For e = 
0.8 (see Fig. 16(a) | and for uj below 0.55, the orbit Of 
is hyperbolic as well as Oi and O^ : The cell is fully 
mixing. Soon after Of turns elliptic at w « 0.55, it 
disappears around uj « 0.57; the remaining orbits are 
hyperbolic until O^ turns elliptic {uj « 1.55). Finally, at 
u! 1.8, the elliptic O^ merges with the hyperbolic Oi 
(as described by the bifurcation (|24p ) to give an elliptic 
Oi : Only partial mixing is achieved. 



For e = 0.4 (see Fig. |16(b)[ ), the situation is slightly 
different : For low uj, the orbit Of (plain line) is present 
and mostly elliptic. Its residue is higher than one only 
in a small range of uj (around uj — 0.61, see the inset of 
Fig. 16(b)), and only this small domain of uj is suitable 
for complete mixing, since Oi and O^ are hyperbolic 
for small uj. Then, for uj « 0.7, the elliptic Of disap- 
pears, and as long as Oi (dashed line) and O^ (dash- 
dotted line) are hyperbolic, the cell is still fully mixing. 
Around w = 1, Oi turns elliptic and so mixing is only 
partial. Moreover, when it turns back to hyperbolicity 
for UJ w 1.55, the orbit O^ soon becomes elliptic : The 
parameter range available for complete mixing is small. 
Finally, O^ merges with the hyperbolic Oi through the 
bifurcation ((24l) , leaving an elliptic Oi : Full mixing can- 
not be achieved any more. 

Figure [T7| summarizes the residue study with the do- 
main of ellipticity/hyperbolicity of these orbits. The do- 
mains of parameters are noted with letters, in agreement 
with the labeling of Fig. [T31 The gray colored domain 
is the most suitable for complete mixing, since Oi and 
O2 are hyperbolic, while Of and Of do not exist. The 
domain (f) would also be suitable; however, in this range 
of parameters, new orbits are born when e decreases. 
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FIG. 15: Modulus of the eigenvalues of orbit Oi for the stream 
function (fTi)) with lu = 1.67. 






FIG. 16: Residue curves for the stream function (|14|l as func- 
tions of the frequency uj : The plain line corresponds to the 
residues of Oi, the dashed one to the residues of and the 
dash-dotted to the ones of (upper curve) and O2 (lower 
curve) . The dotted lines indicate the locations of the bifurca- 
tions. 



IV. CONCLUSION AND PERSPECTIVES 



Time-periodic perturbations are able to generate 
chaotic mixing in two-dimensional channels. Optimal 
mixing is obtained when all regular structures are bro- 
ken by the perturbation. We have shown here that this 
optimal mixing can be obtained in a small domain of 
phase space by combining two strategies : First, by con- 
structing two transport barriers, confining the fluid in a 
bounded region of the channel. Generically, it is expected 
that the motion inside this bounded cell is a mixture of 
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FIG. 17: Mixing domains for the stream function (|14|l in 
parameter space : The letters (a-h) refer to mixing regimes 
as depicted in the typical Poincare sections in Fig. 1131 The 
gray domain represents the expected mixing regime where the 
studied periodic orbits are all hyperbolic (or non-existent). 
The plain curve is associated with the orbit d where its 
residue is equal to one; the dashed curve is associated with 
the orbit O2', the upper dash-dotted line bounds the domain 
of existence of Of (it is elliptic above and does not exist 
below), while the lower dash-dotted line delimits the one of 
of (existing below but not above). While the latter orbit 
is mainly elliptic, the dotted line encloses the small domain 
where it is hyperbolic. 



regular (non-mixing) regions and chaotic (mixing) ones. 
Then, using the linear stability of a few selected periodic 
orbits (represented by their residues) and the identifica- 
tion of bifurcations, we gave conditions on the parameters 
of the system (represented by the amplitude and the fre- 
quency of the periodic forcing) such that a high mixing 
occurs in the cell. We have shown that complete mixing 
is expected in a large region of parameter space. The 
mixing properties have been a posteriori analyzed and 
confirmed by computing a finite time Lyapunov map of 
initial conditions space. Moreover a strategy for placing 
an initial drop of dye has been proposed by combining 
the information of the finite time Lyapunov map with a 
map in initial conditions space of average return times. 
Finally, we have shown that the strategy we developped 
is robust to several effects like truncations of the Fourier 
series giving the exact shape of the perturbation, three- 
dimensional effects and molecular diffusion. When no- 
slip boundary conditions apply, it has to be taken into 
account in the computation of the perturbation since the 
restored transport barrier are not robust without this, 
although the good mixing properties do not seem to be 
affected. 
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